1. Introducción

Una política pública con impacto territorial, ¿tendrá el mismo impacto en el centro de la Ciudad de Buenos Aires, Córdoba o Santa Fe qué en un barrio de una ciudad pequeña? ¿Cómo podemos conocer de antemano que hay en el territorio? Cuando nos enfrentamos a diferentes problemas de políticas públicas, es importante entender y tener en cuenta el territorio donde se van a implementar. Un plan de acción diseñado para la Ciudad de Buenos Aires difícilmente sea replicable en una ciudad del interior debido a diferencias intrínsecas al territorio y la población de uno y otro lugar. Para analizar estas diferencias y estimar el impacto de antemano se puede utilizar las unidades geoestadísticas producidas por INDEC. Una de esas unidades geoestadísticas son los radios censales.

Este documento esta organizado de la siguiente manera. En la Sección 2, vamos a introducir el concepto de Radio Censal, cómo se construye y cómo podemos graficarlos fácilmente. En la Sección 3 vamos a analizar qué información hay disponible para los Radios Censales y vamos a realizar distintas operaciones en conjunto a los radios Censales. En la Sección 4 vamos a utilizar la herramienta GeoRef y normalizar direcciones, nombre de departamentos y, en algunos casos, obtener Latitud y Longitud. Cada sección tiene una parte conceptual y una parte de código. A lo largo del documento, nuestro ejemplo será el Municipio de Lomas de Zamora, en la Provincia de Buenos Aires.

2. Introducción a Radios Censales

2.1. ¿Qué es un Radio Censal y cómo se construye?

Es una unidad geográfica que agrupa entre 100 y 300 viviendas, dependiendo de si son urbanos, rurales o rurales mixtos. En la Argentina, la institución responsable de diseñar y mantener los radios censales es el INDEC. En total hay más de 52.000 radios censales en el país, siendo la Provincia de Buenos Aires el que mas tiene. El código del radio censal viene dado por un numero de código de 9 números, como el siguiente:

\[064900302\]

El código es el resultado de una concatenación de códigos:

  • Los primeros dos dígitos (06) pertenecen a las “Jurisdicciones de primer orden” o División Político-Territorial: las 23 provincias y la Ciudad Autónoma de Buenos Aires. En nuestro ejemplo, se refiere a la Provincia de Buenos Aires.
  • Los tres dígitos siguientes (490) son las “Jurisdicciones de segundo orden” o “Divisiones Político-Administrativas” y son los partidos de la Provincial de Buenos Aires, las comunas de la Ciudad Autónoma de Buenos Aires y los departamentos de las demás provincias. En nuestro ejemplo, corresponde al Municipio de Lomas de Zamora.
  • Los próximos 2 dígitos (03) corresponden a las fracciones urbanas. Es una medida geográfica que agrupa, en promedio, 5000 viviendas.
  • Los últimos dos dígitos (02) corresponden al radio censal.

2.2. ¿Qué información se tiene para los Radios Censales?

En la Argentina hay mucha información públicada por distintos organismos públicos que vienen agregados a nivel de Radios Censales. Por ejemplo, junto a los radios censales, INDEC pública cinco variables básicas de las tres unidades de relevamiento del censo (población, hogares y viviendas): población total por sexo, total de hogares, total de viviendas particulares y total de viviendas particulares habitadas.

A su vez, mucha información social, económica y demográfica proveniente del último censo (2010) que pública el INDEC lo hace con datos agregados en radios censales. Sin embargo, utilizando técnicas de Sistemas de Información Geográfica (SIG o GIS en ingles) y utilizando distintas fuentes, se pueden construir una amplia variedad de capas de información: cantidad de escuelas por radio censal, porcentaje de de personas con necesidades básicas insatisfechas, porcentaje del radio censal ocupado por un Barrio Popular, etc.

2.3. Parte Práctica

Lo primero que vamos a tener que hacer es bajarnos el archivo “shapefile” del territorio con el que vamos a trabajar. La fuente oficial en la Argentina es el INDEC y podemos acceder a ellos a través de la siguiente página. Para este tutorial vamos a trabajar con datos de la Provincia de Buenos Aires, por lo tanto vamos a bajar el siguiente archivo:

https://www.indec.gov.ar/ftp/cuadros/territorio/codgeo/Codgeo_Buenos_Aires_con_datos.zip

Una vez que lo hemos descargado, vamos a abrir el archivo utilizando el paquete sf1.

# Instalamos los paquetes si no están instalados en nuestras computadoras
# install.packages("tidyverse")
# install.packages("sf")

# Carga de los paquetes a utilizar
library(tidyverse)
library(sf)

# Abriendo el shp
radios.buenos.aires <- st_read("Datos/Codgeo_Buenos_Aires_con_datos/", layer = "Buenos_Aires_con_datos")
## Reading layer `Buenos_Aires_con_datos' from data source `/Users/lauti/Google Drive/Lauti/Proyectos/03.Accesibilidad/Datos/Codgeo_Buenos_Aires_con_datos' using driver `ESRI Shapefile'
## Simple feature collection with 19577 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 3721440 ymin: 5452221 xmax: 4335413 ymax: 6305225
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=-90 +lon_0=-66 +k=1 +x_0=3500000 +y_0=0 +ellps=WGS84 +units=m +no_defs

Cuando cargamos el archivo geografico, vemos

  • Simple feature collection with 19577 features and 8 fields: estamos abriendo un dataset que contiene 19577 filas y 8 columnas.
  • geometry type: MULTIPOLYGON: los archivos con información geográfica contienen colecciones de puntos, de líneas, o de polígonos y sus respectivas versiones “múltiples”: múltiples puntos, líneas o polígonos Los Radios Censales son múltiples polígonos
  • dimension: XY: es la cantidad de dimensiones con las que se esta trabajando, en este caso, dos. Si hubiese tres dimensiones, sería XYZ
  • bbox: xmin: 3721440 ymin: 5452221 xmax: 4335413 ymax: 6305225: bbox es la simplificación para “bounding box”, una caja delimitadora. Estos valores son la latitud mínima, la longitud mínima, la latitud máxima y la longitud máxima del conjunto de datos
  • epsg (SRID): NA y proj4string: +proj=tmerc +lat_0=-90 +lon_0=-66 +k=1 +x_0=3500000 +y_0=0 +ellps=WGS84 +units=m +no_defs: nos brindan información referida al sistema de coordenadas de referencia (Coordinate Reference System (CRS)).

Como vemos en el punto anterior, epsg (SRID): es NA. Por lo tanto, lo que se debe hacer es transformarlo al mismo CRS que vamos a utilizar a lo largo del trabajo EPSG: 4326 y para ello vamos a utilizar la función st_transform.

radios.buenos.aires <- st_transform(radios.buenos.aires, crs = 4326)

La primer aproximación a este dataset debería ser ver que datos tiene:

# Exploración inicial de datos
summary(radios.buenos.aires)
##    toponimo_i            link           varon            mujer       
##  Min.   :284416   060070101:    1   Min.   :   0.0   Min.   :   0.0  
##  1st Qu.:299826   060070102:    1   1st Qu.: 196.8   1st Qu.: 217.0  
##  Median :304911   060070103:    1   Median : 383.0   Median : 418.0  
##  Mean   :303920   060070104:    1   Mean   : 388.9   Mean   : 410.2  
##  3rd Qu.:309949   060070105:    1   3rd Qu.: 545.0   3rd Qu.: 577.0  
##  Max.   :336781   060070106:    1   Max.   :3647.0   Max.   :3108.0  
##                   (Other)  :19571   NA's   :25       NA's   :25      
##    totalpobl         hogares       viviendasp       viv_part_h    
##  Min.   :   0.0   Min.   :   0   Min.   :   0.0   Min.   :   0.0  
##  1st Qu.: 417.0   1st Qu.: 148   1st Qu.: 203.0   1st Qu.: 140.0  
##  Median : 801.0   Median : 261   Median : 291.0   Median : 243.0  
##  Mean   : 799.2   Mean   : 245   Mean   : 275.1   Mean   : 226.3  
##  3rd Qu.:1121.0   3rd Qu.: 341   3rd Qu.: 371.0   3rd Qu.: 317.0  
##  Max.   :5895.0   Max.   :1370   Max.   :1259.0   Max.   :1109.0  
##  NA's   :25       NA's   :25     NA's   :25       NA's   :25      
##           geometry    
##  MULTIPOLYGON :19577  
##  epsg:4326    :    0  
##  +proj=long...:    0  
##                       
##                       
##                       
##

Aquí vemos que el archivo radios.buenos.aires tiene las siguientes variables:

  • toponimo_i: Identificación única propia del archivo
  • link: ID del Radio Censal
  • varon: cantidad de varones en dicho radio censal
  • mujer: cantidad de mujeres en dicho radio censal
  • totalpobl: población total
  • hogares: cantidad total de hogares
  • viviendasp: total de viviendas particulares
  • viv_part_h: total de viviendas particulares habitadas
  • geometry: información geográfica que utiliza R para hacer la visualización del mapa.

El siguiente paso será visualizar los radios censales de la Provincia de Buenos Aires. Para dicha tarea hay una gran variedad de paquetes en R que responden a distintas necesidades. Nosotros vamos a utilizar el paquete ggplot2 para las visualizaciones estáticas (preferible, por ejemplo, para la elaboración de informes impresos) y el paquete leaflet para la visualización dinámica (preferible para la exploración interactiva).

Vayamos primero con ggplot2:

# Instalamos los paquetes si no están instalados en nuestras computadoras
# install.packages("ggplot2")

# Carga de los paquetes a utilizar
library(ggplot2)

ggplot() + geom_sf(data = radios.buenos.aires) +
  labs(title = "Radios Censales de la Provincia de Buenos Aires",
         subtitle = "Fuente: INDEC")

Lo primero que se puede ver a simple vista es el contorno de la Provincia de Buenos Aires. Lo segundo que se puede apreciar es el tamaño dispar de muchos de los radios censales, concentrándose los mas pequeños alrededor de las ciudades y los más extensos en el medio del campo. La desventaja es que si como analístas queremos explorar un área con mayor detenimiento, deberíamos usar algún paquete interactivo, como leaflet.

# Instalamos los paquetes si no están instalados en nuestras computadoras
# install.packages("leaflet")

# Carga de los paquetes a utilizar
library(leaflet)

leaflet() %>%
     addTiles() %>%
     addProviderTiles("CartoDB.Positron") %>%
     # Mapeando los poligonos de Lomas de Zamora
     addPolygons(data = radios.buenos.aires,
                 color = "grey", weight = 1, smoothFactor = 0.5,
                 opacity = 1.0)

De esta manera tenemos el mismo mapa que antes, pero ahora podemos explorar el mapa interactivamente, hacer zoom en alguna ciudad o movernos por la provincia.

3. Construyendo Información sobre Radios Censales

En el punto, se habló de que utilizando distintas fuentes y formatos de datos, es posible enriquecer la información que se posee sobre los radios censales. Pero, ¿cuáles son esas fuentes de datos y sus formatos?¿Qué técnicas se pueden utilizar para ellos?

Los SIG modelan la realidad territorial para convertirla en datos geográficos que son manipulados en un entorno informatizado. Para ello utilizan los modelos de representación raster y vectorial.

3.1. Archivos Vectoriales

El modelo de representación vectorial, modeliza las datos valiéndose de primitivas geométricas, tales como puntos, líneas y polígonos. Adosados a dichas geometrías, se encuentran los atributos temáticos de los fenómenos que representan. Por ejemplo, en el caso de los cursos de agua, modelizados a través de polilíneas, se pueden encontrar atributos como el nombre y categoría de los cursos de agua, el régimen hídrico, el caudal anual, etc2.

3.1.1 Archivos planos

Estos archivos normalmente son .xls, .csv, etc. y ya vienen con una columna que hace referencia a la geografía con la que estamos queriendo trabajar. INDEC, por ejemplo, pública información sobre de este la provisión de agua, cloacas, etc. Cuando trabajamos con este tipo de información, para integrarla a nuestros mapas es simplemente cuestión de agregarla.

Vamos a trabajar con los datos del INDEC (en base al Censo 2010) de “Hogares con red pública de agua”. Vamos a descargar los datos directamente desde la web3.

# Definiendo la URL
agua.url <- "https://raw.githubusercontent.com/lautarocantar/Infrastructura_Geografica/master/01_archivos_planos/agua/agua_radiocensal.csv"

# Abriendo el dataset
agua <- read_csv(agua.url,
                 col_types = cols(
                   Link = col_character(),
                   Nombre_de_provincia = col_character(),
                   Total_de_hogares_2010 = col_integer(),
                   Red_pública_de_agua = col_integer(),
                   Perforacion_con_bomba_a_motor = col_integer(),
                   Perforacion_con_bomba_manual = col_integer(),
                   Pozo = col_integer(),
                   Transporte_por_cisterna = col_integer(),
                   Agua_delluvia_rio_o_canal = col_integer(),
                   Total = col_integer(),
                   Total_sin_red = col_integer(),
                   Porcentaje_dehogares_con_red_de_agua = col_double(),
                   Porcentaje_de_hogares_sin_red_de_agua = col_double())
                 )

Nuevamente, lo que debemos hacer antes de empezar a utilizar los datos para un mapa, debemos analizar qué variables tenemos, cuáles son los valores minimos, maximos, etc. Como vemos con los resultados de summary(agua), hay 52406 radios censales relevados y 11 variables.

summary(agua)
##  Radio_Censal        Provincia         Cantidad_hogares
##  Length:52406       Length:52406       Min.   :   0.0  
##  Class :character   Class :character   1st Qu.: 134.0  
##  Mode  :character   Mode  :character   Median : 245.0  
##                                        Mean   : 232.2  
##                                        3rd Qu.: 322.0  
##                                        Max.   :1517.0  
##                                                        
##  agua_conexion_red_publica agua_perforacion_con_bomba_a_motor
##  Min.   :   0.0            Min.   :   0.00                   
##  1st Qu.:  26.0            1st Qu.:   0.00                   
##  Median : 217.0            Median :   1.00                   
##  Mean   : 194.9            Mean   :  27.06                   
##  3rd Qu.: 301.0            3rd Qu.:  12.00                   
##  Max.   :1315.0            Max.   :1176.00                   
##                                                              
##  agua_perforacion_con_bomba_manual agua_conexion_a_pozo
##  Min.   :  0.000                   Min.   :  0.000     
##  1st Qu.:  0.000                   1st Qu.:  0.000     
##  Median :  0.000                   Median :  0.000     
##  Mean   :  1.635                   Mean   :  4.867     
##  3rd Qu.:  1.000                   3rd Qu.:  3.000     
##  Max.   :199.000                   Max.   :623.000     
##                                                        
##  agua_transporte_por_cisterna agua_agua_de_lluvia_rio_o_canal
##  Min.   :  0.000              Min.   :  0.000                
##  1st Qu.:  0.000              1st Qu.:  0.000                
##  Median :  0.000              Median :  0.000                
##  Mean   :  1.736              Mean   :  2.031                
##  3rd Qu.:  0.000              3rd Qu.:  0.000                
##  Max.   :373.000              Max.   :374.000                
##                                                              
##  agua_total_sin_conexion_a_red agua_pje_hogares_red
##  Min.   :   0.00               Min.   :0.0000      
##  1st Qu.:   0.00               1st Qu.:0.3137      
##  Median :   3.00               Median :0.9856      
##  Mean   :  37.33               Mean   :0.7139      
##  3rd Qu.:  28.00               3rd Qu.:1.0000      
##  Max.   :1224.00               Max.   :1.0000      
##                                NA's   :128

Las variables que tenemos son:

  • Radio_Censal: Radio Censal
  • Provincia: Nombre de la provincia
  • Cantidad_hogares: Cantidad de hogares totales en dicho Radio Censal
  • agua_conexion_red_pública: Cantidad de hogares con conexión a red pública
  • agua_perforacion_con_bomba_a_motor: Cantidad de hogares con perforaciones con bomba a motor
  • agua_perforacion_con_bomba_manual: Cantidad de hogares con perforaciones con bomba a manual
  • agua_conexion_a_pozo: Cantidad de hogares con conexión a pozo
  • agua_transporte_por_cisterna: Cantidad de hogares con provisión de agua por transporte por cisterna
  • agua_agua_de_lluvia_rio_o_canal: Cantidad de hogares con agua de lluvia, rio o canal
  • agua_total_sin_conexion_a_red: Cantidad de hogares sin conexión a red (sumatoria de todas las variables anteriores, salvo “agua_conexion_red_pública”
  • agua_pje_hogares_red: Porcentaje de los hogares de cada Radio Censal que tienen conexión a la red de agua pública.

Antes de hacer la unión entre el dataset agua y el dataset de radios.buenos.aires, vamos a hacer una subselección para trabajar únicamente con los datos del Municipio de Lomas de Zamora. Para eso vamos a seleccionar aquellos radios censales que comiencen con los valores “06490”. En Lomas de Zamora hay 612 Radios Censales.

# Subseleccion de # Subseleccion de Radios Censales en Lomas de Zamora
agua.lomas.zamora <- agua %>% filter(str_detect(Radio_Censal, "06490"))

# Subseleccion de Radios Censales en Lomas de Zamora
rc.lomas.zamora <- radios.buenos.aires %>% filter(str_detect(link, "06490"))

Ahora si, podemos hacer la unión de ambas bases de datos:

# Cambiando el tipo de la variable "Link" para hacer la union de los datos
radios.buenos.aires$link <- as.character(radios.buenos.aires$link)

# Union de ambas bases de datos
rc.lomas.zamora <- rc.lomas.zamora %>% 
  left_join(agua %>% select(Radio_Censal, agua_pje_hogares_red), by = c("link" = "Radio_Censal"))

Perfecto, ya tenemos la unión de ambos dataset y debemos graficarlo. Pero antes vamos a explorar cómo están distribuidos los datos del porcentaje de los hogares de cada Radio Censal que tienen conexión a la red de agua pública. Para eso vemos los cuantiles de agua_pje_hogares_red.

quantile(rc.lomas.zamora$agua_pje_hogares_red)
##        0%       25%       50%       75%      100% 
## 0.0000000 0.9787033 0.9919514 0.9974109 1.0000000

Dado que Lomas de Zamora es enteramente urbano, tiene un altiísimo porcentaje de hogares con conexión a la red puública. Por lo tanto, para graficarlo vamos a generar una nueva variable que clasifique a los radios censales en 4 categorías diferentes, una categoriía para cada cuartil.

# Nueva Variable
rc.lomas.zamora <- rc.lomas.zamora %>% 
  mutate(agua = case_when(.$agua_pje_hogares_red <= 0.9787033 ~ 1,
                          .$agua_pje_hogares_red <= 0.9919514 ~ 2,
                          .$agua_pje_hogares_red <= 0.9974109 ~ 3,
                          .$agua_pje_hogares_red > 0.9974109 ~ 4)
         )

Vamos a crear una paleta de colores para cada una de las categorías para poder graficarlas y que se entienda nuestro hermoso gráfico.

# Cambiando los nombres de las variables creadas anteriormente
cuartiles <- factor(c(1,2,3,4), 
                    labels = c("Primer cuartil", "Segundo cuartil", "Tercer cuartil", "Cuarto cuartil"))

# Creando la paleta de colores
Color_Assets <- colorFactor(c("#e3f2fd", "#64b5f6", "#1e88e5", "#0d47a1"),
                           levels = cuartiles, ordered=TRUE)

binpal <- colorBin(colors, rc.lomas.zamora$agua, 4, pretty = FALSE)

Ahora si, a hacer el mapa!

# Mapeando
leaflet() %>%
  addTiles() %>%
  addProviderTiles("CartoDB.Positron") %>%
  # Mapeando los poligonos de Lomas de Zamora
  addPolygons(data = rc.lomas.zamora,
              color = ~binpal(agua), weight = 1, smoothFactor = 1,
              stroke = FALSE, fillOpacity = 0.5) %>% 
  addLegend("bottomright", pal = Color_Assets, values =  cuartiles,
    title = "Porcentaje de Acceso <br>a la Red Pública",
    opacity = 1
  )

3.1.2. Archivos geográficos de punto

Cada lugar, cada movimiento que hacemos en la calle, tiene una latitud y una longitud asociada. Esa latitud y longitud es lo que después se utiliza en el mapa para localizar la ubicación en el lugar pertinente. Cuando tenemos esa latitud y longitud, podemos hacer varias operaciones para agregarlas a Radios Censales. Una de esas operaciones es la intersección entre los polígonos de los radios censales y los puntos que estamos interesados. De esta forma, todos los puntos que caigan dentro de dicho polígono serán asignados a dicho radio censal.

Para nuestro ejemplo vamos a trabajar con un dataset de las Farmacias que atienden a la obra social PAMI en el Municipio de Lomas de Zamora. Los vamos a abrir desde la URL y vemos la distribucion de los primeras observaciones4.

# Definiendo la URL
farmacias.url <- "https://raw.githubusercontent.com/lautarocantar/Infrastructura_Geografica/master/pami.lanus.csv"

# Abriendo el dataset
farmacias <- read_csv(farmacias.url, 
                      col_names = c("Provincia", "UGL_ID", "UGL_DESC", "Descripcion", "Direccion", "Localidad",
                                    "Codigo_Postal", "Longitud", "Latitud"), skip = 1)

# Viendo los datos
head(farmacias, 10) %>% knitr::kable()
Provincia UGL_ID UGL_DESC Descripcion Direccion Localidad Codigo_Postal Longitud Latitud
BUENOS AIRES 10 LANUS GANDULFO BALCARCE 00314 LANUS B1832BFT -58,392801245 -34,762076531
BUENOS AIRES 10 LANUS ARTAZA JACARANDA 00307 LANUS NA -58,379397510 -34,689798551
BUENOS AIRES 10 LANUS CAMPI MAROTTO 00279 LANUS B1806DYP -58,568127673 -34,893897918
BUENOS AIRES 10 LANUS AMERICANA GARCIA 00202 LANUS NA 0,000000000 0,000000000
BUENOS AIRES 10 LANUS FERNANDEZ AV MITRE 02044 LANUS B1870ESA -58,353043000 -34,670996000
BUENOS AIRES 10 LANUS CAVA MARCONI 00798 CLAYPOLE B1849FEA -58,347344367 -34,806226408
BUENOS AIRES 10 LANUS DEL PINO FONROUGE Y AV ALS 00000 LANUS NA 0,000000000 0,000000000
BUENOS AIRES 10 LANUS FARMA EXPRESS SCS BV GOB MARTIN RODR 01850 LANUS NA -58,373612041 -34,705055918
BUENOS AIRES 10 LANUS FERNANDEZ AV GALICIA 00054 LANUS B1870EVM -58,382771735 -34,677032061
BUENOS AIRES 10 LANUS PROTTO EL CEIBO 00199 LANUS B1842ECD -58,496714000 -34,832442000

La primer tarea que tenemos es conventir el archivo .csv en un archivo espacial. Para eso vamos a utilizar el paquete rgdal que nos permite hacer algunas operaciones sobre archivos geográficos. Luego creamos una copia del archivo farmacias para alterar radicalmente el archivo (lo vamos a usar en la próxima sección).

# Instalamos los paquetes si no están instalados en nuestras computadoras
# install.packages("rgdal")

# Carga de los paquetes a utilizar
library(rgdal)

# Copia del dataset
farmacias.shp <- farmacias

Si quisieramos hacer la transformación con los datos como están, nos dariía error, daddo que los datos no estaán en formato numeérico. Por lo tanto, para hacer esa operación, primero debemos reemplazar las comas por los puntos en las columnas Latitud y Longitud y luego si convertiras a formato numérico.

# Reemplanzando las comas por los puntos en las columnas Latitud y Longitud
farmacias.shp$Latitud <- gsub(",", '.', farmacias.shp$Latitud, fixed = T)
farmacias.shp$Longitud <- gsub(",", '.', farmacias.shp$Longitud, fixed = T)

# Convirtiendolas en formato numérico
farmacias.shp$Latitud <- as.numeric(farmacias.shp$Latitud)
farmacias.shp$Longitud <- as.numeric(farmacias.shp$Longitud)

Ahora si estamos en condiciones para convertir el archivo desde un .csv a un SpatialPointsDataFrame, utlilizando la función coordinates.

# Transformando el formato del archivo
coordinates(farmacias.shp) <- ~Longitud+Latitud
farmacias.shp <- st_as_sf(farmacias.shp, crs = 4326, coords = c("Longitud", "Latitud"))

Veamos a simple vista como quedarían mapeadas dichas farmacias en los Radios Censales del Departamento de Lomas de Zamora. Tal como se ve en el mapa debajo, hay muchas farmacias ubicadas por fuera del territorio de Lomas de Zamora (inclusive por fuera del territorio nacional, los cuales son error de carga de los datos).

leaflet() %>%
  addTiles() %>%
  addProviderTiles("CartoDB.Positron") %>%
  # Mapeando los polígonos de Lomas de Zamora
  addPolygons(data = rc.lomas.zamora,
              color = "grey", weight = 1, smoothFactor = 1,
              stroke = 0.1, fillOpacity = 0.5) %>%
  # Mapeando los puntos de las Farmacias
  addCircleMarkers(data = farmacias.shp, 
    radius = 2,
    color = "red",
    stroke = FALSE, fillOpacity = 0.5)

Antes de hacer la interseccioón, debemos considerar que ambos archivos deben tener la misma proyección:

# Controlando que tengan la misma proyeccion
st_crs(rc.lomas.zamora)
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(farmacias.shp)
## Coordinate Reference System: NA

Como nuestro archivo farmacias.shp no tiene asignada una CRS, se la asignamos

st_crs(farmacias.shp) = 4326

Y ahora si podemos hacer la intersección

# Creando la interseccion
inter.farmacias <- st_intersection(rc.lomas.zamora, farmacias.shp)
inter.farmacias <- inter.farmacias %>% group_by(link) %>% tally() %>% as.data.frame()

rc.lomas.zamora <- rc.lomas.zamora %>% left_join(inter.farmacias %>% select(link, n))

pal <- colorNumeric(
  palette = "Blues",
  domain = rc.lomas.zamora$n)


leaflet() %>%
  addTiles() %>%
  addProviderTiles("CartoDB.Positron") %>%
  # Mapeando los polígonos de Lomas de Zamora
  addPolygons(data = rc.lomas.zamora,
              color = ~pal(n), weight = 1, smoothFactor = 1,
              stroke = 0.1, fillOpacity = 0.5)

3.1.3. Archivos de polígonos

Otra fuente de información geográfica que podemos utilizar son otros polígonos construidos para otros objetivos como, por ejemplo, la zonificación de las áreas de una ciudad o un relevamiento de barrios populares. Por ejemplo, un radio censal puede estar atravesado por uno o más barrios populares o tener una o más clasificaciones zonales. Una estrategia para trabajar con este tipo de situaciones es utilizar nuevamente la operación intersección y luego calcular el porcentaje de superposición entre ambos polígonos.

Para este ejemplo, vamos a seguir trabajando con el Departamento de Lomas de Zamora y vamos a sumar el dataset de Barrios Populares, un relevamiento en territorio realizado por distintas organizaciones sociales y Jefatura de Gabinete de Ministros a diciembre 2016. El mapa nos muestra poliígonos de los los barrios populares de Argentina. La información la podés descargar desde este link.

# Abriendo el shp
barrios.populares <- st_read(dsn="Datos/barrios-populares/", layer = "barrios-populares")
## Reading layer `barrios-populares' from data source `/Users/lauti/Google Drive/Lauti/Proyectos/03.Accesibilidad/Datos/barrios-populares' using driver `ESRI Shapefile'
## Simple feature collection with 4228 features and 5 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -71.56234 ymin: -54.82007 xmax: -54.01357 ymax: -22.04006
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
# Filtrando para Buenos Aires
barrios.populares <- barrios.populares %>% filter(departamen == "LOMAS DE ZAMORA")

Haciendo una una vista muy rápida de los radios censales y los barrios populares que ubicados en Lomas de Zamora vemos que claramente hay una superposición entre ambos polígonos.

leaflet() %>%
  addTiles() %>%
  addProviderTiles("CartoDB.Positron") %>%
  # Mapeando los polígonos de Lomas de Zamora
  addPolygons(data = rc.lomas.zamora,
              color = "grey", weight = 1, smoothFactor = 1,
              stroke = 0.1, fillOpacity = 0.5) %>%
  # Mapeando los polígonos de los Barrios Populares
    addPolygons(data = barrios.populares,
              color = "blue", weight = 1, smoothFactor = 1,
              stroke = 0.1, fillOpacity = 0.5)

Una forma de integrar ambas bases de datos es calcular que porcentaje de un radio censal está ocupada por un barrio popular. Para esto, primero tenemos que calcular el área de cada radio censal y de cada barrio popular.

# Calculando el area en KM2 de cada unos de los poligonos
rc.lomas.zamora$Area_RC <- st_area(rc.lomas.zamora)
barrios.populares$Area_BP <- st_area(barrios.populares)

Antes de hacer la intersección entre ambos archivos, nos tenemos que asegurar que tengan la misma proyección.

# Controlando que tengan la misma proyeccion
st_crs(rc.lomas.zamora)
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(barrios.populares)
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

Ahora vamos a hacer intersección entre ambos archivos. Para eso vamos a utilizar la función st_intersection y crear un nuevo archivo que se llame int. Luego vamos a calcular el área (en \(m^2\)) de la superposicioón.

# Creando la interseccion
int <- st_intersection(rc.lomas.zamora, barrios.populares)

# Calculando el area de la superposicion 
int$area_int <- st_area(int)

El paso siguiente es calcular el porcentaje de la superpocisioón del aárea, la cual no puede ser mayor a 100.

# Calculando el porcentaje del RC que es un BP
int$Pje_BP_en_RC <- (int$area_int/int$Area_RC)*100

# Explorando los resultados
summary(int$Pje_BP_en_RC)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##   0.00002   0.55638  11.29729  30.54083  60.26615 100.00000

Por uúltimo hacer un par de operaciones para acomodar los datos y adjuntarnos a nuestro dataset rc.lomas.zamora

# Operaciones intermedias 
int.2 <- as.data.frame(int) %>% group_by(link) %>% summarise(Pje_Total_BP_en_RC = sum(Pje_BP_en_RC))
int.2$Pje_Total_BP_en_RC <- round(int.2$Pje_Total_BP_en_RC, 2)

# Adjuntarlo al nuestro dataset original
rc.lomas.zamora <- rc.lomas.zamora %>% left_join(int.2)
rc.lomas.zamora$Pje_Total_BP_en_RC <- as.numeric(rc.lomas.zamora$Pje_Total_BP_en_RC)

Por último graficamos aquellos radios censales que tengan Barrios Populares:

pal <- colorNumeric(
  palette = "Blues",
  domain = rc.lomas.zamora$Pje_Total_BP_en_RC)


leaflet() %>%
  addTiles() %>%
  addProviderTiles("CartoDB.Positron") %>%
  # Mapeando los polígonos de Lomas de Zamora
  addPolygons(data = rc.lomas.zamora,
              color = ~pal(Pje_Total_BP_en_RC), weight = 1, smoothFactor = 1,
              stroke = FALSE, fillOpacity = 0.5) %>%
  addLegend("bottomright", pal = pal, values =  rc.lomas.zamora$Pje_Total_BP_en_RC,
    title = "Porcentaje de Barrio Popular<br>en Radios Censales",
    opacity = 1
  )

3.2. Archivos raster

Los archivos raster constan de una matriz de celdas (o píxeles) organizadas en filas y columnas (o una cuadrícula) en la que cada celda contiene un valor que representa información, como la temperatura. Los rásteres son fotografías aéreas digitales, imágenes de satélite, imágenes digitales o incluso mapas escaneados5. Para poder llevar esta información a nuestros radios censales, se requieren operaciones más complejas, por eso no lo vamos a ver en este momento.

Una aplicacion muy interesante de este tipo de archivos geograficos es el trabajo que llevan adelante desde la empresa Dymaxion Labs con la detección de Asentamientos Precarios. Utilizando imagenes satelitales y algorítmos de aprendizaje automático, identifican patrones en los techos de las casas y distinguen que zonas presentan patrones similares a los de asentamientos precarios.

4. GeoRef

4.1. Conociendo a GeoRef

Arranquemos esta sección haciendo un breve ejercicio: agarremos un lápiz y papel y durante un minuto escribamos todas las variantes de la localidad, departamento o provincia que elijamos. Esta es mi lista:

  • Buenos Aires
  • Prov. de Bs. As.
  • PBA
  • Bs. As.
  • B. Aires

Para el ojo humano, es obvio que todas las variantes hacen referencia a la Provincia de Buenos Aires pero cuando trabajamos con sistemas informáticos, son todas provincias diferentes. Esto genera un nuevo desafío: ¿cómo hacemos para que nuestro sistema informático lo reconozca como lo mismo? Para solucionar este problema, desde datos.gob.ar se creo una herramienta que permite normalizar y codificar los nombres de unidades territoriales de la Argentina (provincias, departamentos, municipios y localidades) y de sus calles, así como ubicar coordenadas dentro de ellas.

El Servicio de Normalización de Datos Geográficos es una API REST. Una API es una interfaz de programación de aplicaciones (la sigla API viene del ingles Application Programming Interface). Es un conjunto de rutinas que provee acceso a funciónes de un determinado software. Supongamos que queremos consultar los departamentos de la provincia de Córdoba. A través de una consulta web como la que se ve debajo, la API nos devuelve cierta información. Primero desglocemos la URL que estamos utilizando6. Notemos que los parametros empiezan despues del ? y se separan por &.

\[https://apis.datos.gob.ar/georef/api/departamentos?formato=csv&provincia=cordoba&max=10\]

  • https://apis.datos.gob.ar/georef/api/: dirección básica
  • departamentos: unidad geográfica que nos intereasa
  • ?formato=csv: formato en qué queremos que nos brinde los dato
  • &provincia=cordoba: provincia de la que estamos interesados consultar
  • &max=10: cantidad máxima de registro que nos interesa (esto se puede cambiar)

El resultado es el siguiente:

departamento_id departamento_nombre departamento_centroide_lat departamento_centroide_lon provincia_id provincia_nombre
14161 Tercero Arriba -32.28771 -63.77925 14 Córdoba
14042 General San Martín -32.51433 -63.25622 14 Córdoba
14049 Ischilín -30.39972 -64.60949 14 Córdoba
14077 Pocho -31.46344 -65.43848 14 Córdoba
14091 Punilla -31.22256 -64.58615 14 Córdoba
14154 Sobremonte -29.76546 -64.14289 14 Córdoba
14021 Colón -31.14405 -64.15287 14 Córdoba
14028 Cruz del Eje -30.65806 -65.07747 14 Córdoba
14056 Juárez Celman -33.32945 -63.60634 14 Córdoba
14098 Río Cuarto -33.33080 -64.49418 14 Córdoba

4.2. Parte Práctica

4.2.1 GeoRefAr para normalizar nombres

Para trabajar con la API de GeoRef, vamos a utilizar el paquete georefar7, desarrollado por Patricio Del Boca. Dado que este paquete no se encuentra disponible en CRAN, hay que utlizar un paquete intermedio llamado devtools (si no esta instalado en tu computadora, es posible que lo tengas que instalar).

# Instalacion del paquete "georefar"
# install.packages("devtools")
# devtools::install_github("pdelboca/georefar")

# Cargamos el paqute "georefar"
library(georefar)

Para este ejemplo vamos a utilizar nuevamente el dataset farmacias. Recordemos que este dataset muestra todas las farmacias de la zona de la Agencia de Lanús. En el ejemplo anterior utilizamos las Latitudes y Longitudes, pero qué pasaría si no tuvieramos dicha información y quisieramos identificar únicamente las farmacias que están localizadas en el Departamento de Lomas de Zamora (o, propiamente dicho, Municipio de Lomas de Zamora). Una de las formas que tenemos de hacer esto es utilizar la variable Localidad, que nos informa en qué localidad están las farmacias. Sin embargo, un departamento posee varias localidades.

Para solucionar esto, podemos utilizar la función get_localidades para Lomas de Zamora y el resultado son las siguientes 7 localidades:

# LLamado a la API utilizando el paquete "georefar"
localidades <- get_localidades(provincia = "Buenos Aires", departamento = "LOMAS DE ZAMORA") 

# Transformacion de los resultados
nombre.localidades <- localidades %>% select(nombre) %>% rename(Loalidad = nombre)
nombre.localidades %>% knitr::kable()
Loalidad
LOMAS DE ZAMORA
LLAVALLOL
BANFIELD
VILLA CENTENARIO
TEMPERLEY
TURDERA
VILLA FIORITO

Por lo tanto, una vez que tenemos los nombres normalizados de las localidades de Lomas de Zamora y las localidades de las farmacias, el proximo paso es simplemente asignarles el nombre del departamento a aquellas farmacias pertenecientes a alguna de las 7 localidades anteriores. Para eso creamos una variable para los casos en los que la localidad sea de Lomas de Zamora.

# Creamos un vector con los nombres de las localidades
nombre.localidades <- as.character(unlist(nombre.localidades[,1]))

# Creamos una variable para los casos en los que la localidad sea de Lomas de Zamora 
farmacias <- farmacias %>% mutate(Departamento = case_when(Localidad %in% nombre.localidades ~ "Lomas de Zamora"))

# Filtramos las farmacias de Lomas de Zamora
farmacias.lomas <- farmacias %>% filter(Departamento == "Lomas de Zamora")
head(farmacias.lomas) %>% knitr::kable()
Provincia UGL_ID UGL_DESC Descripcion Direccion Localidad Codigo_Postal Longitud Latitud Departamento
BUENOS AIRES 10 LANUS CATARINA FARMA SCS GIACHINO 03623 BANFIELD B1828EQS -58,477408616 -34,739879091 Lomas de Zamora
BUENOS AIRES 10 LANUS PFM ITATI SCS VIRGEN DE ITATI 00802 BANFIELD B1828BUZ -58,433920000 -34,722593000 Lomas de Zamora
BUENOS AIRES 10 LANUS MILANO DOMINGO FRENCH 00690 BANFIELD B1828FXE -58,403898653 -34,742336102 Lomas de Zamora
BUENOS AIRES 10 LANUS ALLEGRINO AV MEEKS 01277 TEMPERLEY B1834CGW -58,396823333 -34,776324667 Lomas de Zamora
BUENOS AIRES 10 LANUS LAPINE AV EVA PERON 03305 TEMPERLEY B1834CAQ -58,358549367 -34,760730857 Lomas de Zamora
BUENOS AIRES 10 LANUS NADIR SALUD FARMACIAS ALSINA 00003 BANFIELD B1828KFW -58,391132224 -34,736740612 Lomas de Zamora

Una vez que se llega a este punto, llevar esta informacion al mapa es simplmente cuestion de aplicar una de la técnicas analizadas anteriormente.

4.2.2 GeoRefAr para normalizar direcciones

La API de GeoRef tambien se puede utilizar para normalizar direcciones. Por ejemplo, aquí vamos a tratar de normalizar la dirección “Boedo 202, Lomas de Zamora”.

boedo <- normalizar_direccion("BOEDO 00202, LOMAS DE ZAMORA", departamento = "06490")
boedo %>% knitr::kable()
altura departamento_id departamento_nombre id nombre nomenclatura provincia_id provincia_nombre tipo ubicacion_lat ubicacion_lon
202 06490 Lomas de Zamora 0649001004335 MARIANO BOEDO MARIANO BOEDO 202, Lomas de Zamora, Buenos Aires 06 Buenos Aires CALLE -34.75523 -58.43091
202 06490 Lomas de Zamora 0649001002855 FELIPE BOERO FELIPE BOERO 202, Lomas de Zamora, Buenos Aires 06 Buenos Aires CALLE NA NA

La API funciona de la siguiente manera: recibe un texto y lo divide segun nombre y número. Luego compara el texto con el listado de textos que están en la base de datos y busca mejor aproximación8: en este caso, las majores aproximaciones son BOEDO y BOERO. Luego, busca el número y asigna el segmento de las cuadras. Si tiene Latitud y Longitud, tambien lo provee.

5. Referencias y material extra de consulta

Aquí hay una lista de recursos adicionales para explorar mas de estas cuestiones.

Podes encontrar mas información y datos el repositorio de datos geográficos. Tambien podés explorar el catálogo de datos abiertos en datos.gob.ar y explorar diferentes opciones. Si tenés alguna inquietud o pregunta, no dudes en escribirme a lautaro.cantar.ar@gmail.com. Si te gusto lo que leíste, te invito a compartilo.

  1. Aquí un manual y un tutorial, ambos en inglés↩

  2. Link: http://www.ign.gob.ar/sig↩

  3. En el Repositorio de Datos Geográficos vas a encontrar el código completo de las operaciones sobre el dataset descargado directamente de INDEC↩

  4. El archivo que vamos a utilizar es un subset del archivo original, el cual está disponible acá↩

  5. http://desktop.arcgis.com/es/arcmap/10.3/manage-data/raster-and-images/what-is-raster-data.htm↩

  6. En este post, Natalia Sampietro, Directora de Datos Públicos de la Secretaria de Gobierno de Modernizacion de Argentina, cuenta cosas muy interesantes sobre como funcióna la herramienta. Aquí podés encontrar la documentación sobre la API↩

  7. Copyright (c) 2018 Patricio Del Boca↩

  8. Si la palabra tiene 3 o menos caracteres, la API busca la el nombre exacto, si la palabra tiene entre 4 y 7 caracteres, la API permite una transformación y si tiene más de 8, permite 2 transformaciones.↩